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I. INTRODUCTION 


Previous work at NASA LeRC has shown that flow distortions in aircraft engine inlet ducts can 
be significantly reduced by mounting vortex generators, or small wing sections, on the inside 
surface of the engine inlet f 1], The placement of the vortex generators is an important factor in 
obtaining the optimal effect over a wide operating envelope. In this regard, the only alternative 
to a long and expensive test program which would search out this optimal configuration is a good 
prediction procedure which could narrow the field of search. Such a procedure has been 
developed in collaboration with NASA LeRC, and results obtained by NASA personnel indicate 
that it shows considerable promise for predicting the viscous turbulent flow in engine inlet ducts 
in the presence of vortex generators [1]. 

The prediction tool is a computer code which numerically solves the reduced Navier-Stokes 
equations and so is commonly referred to as RNS3D. Obvious deficiencies in RNS3D have been 
addressed in previous work [3]. Primarily, it is known that the predictions of the mean velocity 
field of a turbulent boundary layer flow approaching separation are not in good agreement with 
data [1]. It was suggested that the use of an algebraic mixing-length turbulence model in RNS3D 
is at least partly to blame for this. Additionally, the current turbulence model includes an 
assumption of isotropy which will ultimately fail to capture turbulence-driven secondary flow 
known to exist in noncircular ducts [30], 

Because turbulent flows contain time and length scales which can vary dramatically, it is 
widely acknowledged that two-equation turbulence models represent the minimum acceptable 
level of closure for general modeling [31]. Such a model must also have the capability to 
integrate the solution all the way to the wall, since flows containing vortex generators will be 
highly dependent on the structure of the turbulence in the near-wall region. Recent work has also 
demonstrated that an anisotropic eddy viscosity, in conjunction with a two-equation turbulence 
model, will significantly enhance the prediction of anisotropic effects [30]. 

For these reasons, the current effort was undertaken to incorporate a low-Reynolds number, 
anisotropic k-e model directly into the RNS3D code, so that the improved modeling capability 
demonstrated in previous work would be accessible for internal flows without sacrificing any of 
the efficiency which RNS3D makes available to its users. 
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II. LOW-REYNOLDS NUMBER TURBULENCE MODELING 


The Governing Equations 


In this work, the two-equation turbulence model based on the turbulent kinetic energy, k, and the 
turbulent energy dissipation rate, e, was developed, implemented, and tested within the context 
of the reduced Navier-Stokes equations as they appear in the RNS3D code. While many 
variations of two-equation turbulence models exist, the k-e formulation has the distinct advantage 
of having been widely used and studied over the longest period of time. 

As with the algebraic mixing-length model, the two-equation model provides for an eddy 
viscosity which is based on characteristic length and velocity scales of the turbulence. In this 
case, the velocity scale is determined by the square root of the turbulent kinetic energy, for which 
a transport equation can be derived exactly from the Navier-Stokes equations using Reynolds- 
averaging. The resulting equation, once some modeling is completed, can be written 
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Here, the time rate of change and convection terms, which appear on the left-hand side of the 
equation, are captured exactly, as are the molecular diffusion, production and dissipation terms 
which appear on the right. Turbulent diffusion effects, which appear as a factor of the turbulent 
viscosity, p„ are approximated by a gradient-diffusion assumption from which the empirical 
coefficient o k arises. Finally, the term represented by D arises from near-wall modeling and 
typically takes the value of the dissipation at the wall [14]. 

The characteristic length scale of the turbulence can be derived from the ratio 

l = k™ / £ (2) 


although a transport equation for e is comparatively more difficult to obtain. Exact equations 
can be derived from the momentum disturbance equations or for enstrophy, the averaged 
magnitude of the fluctuating vorticity, which approaches the dissipation rate for high Reynolds 
numbers. A practical model can then be deduced through order-of-magnitude and dimensional 
arguments. The final form is fairly well-established, and can be written 
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Three empirical coefficients appear in the equation for the turbulent energy dissipation, and they 
are represented by C el , C e2 , and o £ . The term E arises from near-wall modeling, but is typically 
ad-hoc in nature [14]. 

Usually, the Reynolds stresses are represented by the Boussinesq approximation, 
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where the mean strain rate tensor, E {] , is defined as 
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The eddy viscosity, in the k-e formulation, is determined by 
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in which C M is an empirical constant. 

Under the RNS3D assumptions, the streamwise diffusion terms were neglected from both 
equations [4], Incorporating the RNS3D coordinate system and linearizing the equations then 
gives 
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for the ^-equation, where a hat denotes an implicitly treated variable. The production term, 
which depends on 0, the strain-rate invariant, and the low-Reynolds number function D were 
both treated explicitly. Similarly, the e-equation is given by 


h dx 

3 





( \ 

■ 

d£ d& 

+ n v -t- nvu . - 

1 1 d 

U 

ii + 

d€ 

K *dx K °dx 

Re h dx 

i l 

H* 

a 

dx 

1 2 

l 


v e ) 

l 


1 1 d 


Re h dx 



( ^ 

■ 


M' ' 

t i 

dt 

h 

[i + — 

— 


a 

dx 


l e ) 

2 


+ pC C $£ + pC , 

r £l |1 r e2| 


< e' 1 

2 

" € 1 

k 

\ ) 

k - 2pC 
r €2 

Jj 


( 8 ) 


where, once again, the hatted variables denote those which are treated implicitly. The 
linearization of the destruction term here is borrowed from that used in previous studies, where 
it was found to be successful. The low-Reynolds number function E is treated explicitly. 
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The Low-Reynolds Number Modifications 

Despite success with the Launder-Sharma modifications [18] in previous studies [3], it was found 
that the low-Reynolds number function E suggested by Launder-Sharma, when modified 
according to the assumptions in RNS3D, became unstable. For this reason, the low-Reynolds 
number modifications of Chien were adopted [5], These modifications consisted of the functions 

D = 2vk I y 2 , (9) 

which models the wall value of e, allowing the boundary condition on e to be zero, and 

E = (2ve / y 2 )exp(-0.5y *) , (10) 

which also models a '’wall" dissipation effect. Damping functions multiply the empirical 
constants C M and C e2 . They are 

/ = 1 - exp(-0.0115/) (11) 

for the damping of the eddy viscosity, and 

f 2 = 1 - (2/9)exp[-(R r /6) 2 ] , (12) 

which decreases the destruction of dissipation in the near-wall region. Note that, because the 
functions depend on both the distance to the wall, y, and the normalized distance to the wall, y + , 
the normal distance to the wall, as well as the friction velocity at the wall, must be calculated. 
Five empirical constants appear in the model. They are set, according to Chien [5], as 

C = 0.09, C = 1.35, C = 1.8, o = 1.0, o - 1.3. (13) 

ji ’el ’ e2 k e 


Initial Values 

Equations (5)-( 1 3) describe a coupled system of equations which can be solved within the 
RNS3D code by streamwise marching. Thus, in addition to boundary conditions, for example, 
at the wall, a complete description of the turbulent quantities are needed at the inlet of the duct. 
In this study, the values obtained by the existing mixing-length model at the first streamwise 
station are utilized to obtain initial profiles for k and e. The approach taken is equivalent to that 
used in deriving wall-functions for the two variables. That is, it is assumed that production is 
equal to dissipation throughout the inlet plane, allowing k and e to be derived from the mixing 
length eddy viscosity. The resulting equation for k is 

k = C 1/2 / 2 $ (14) 

ti 

where / is the turbulent mixing length and $ is the strain-rate invariant, and 

6 = CV I l. (15) 

n 

Previous experience has shown the utility of damping C M in accordance with the low-Reynolds 
number modifications. In this work, the damping on the initial plane was based on the function 
of Launder-Sharma [18], 
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= exp[-3.4/(l +/? r / 50) 2 ] 

where R T , the turbulent Reynolds number, is based on the mixing length eddy viscosity, 
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When Eq. (17) is substituted into Eq. (16), the implicit relation 

= 0.09exp{-3.4/[l+p/y$/(50pCr)] 2 } ( 18 ) 

is obtained. In this work, Eq. (18) was solved using no more than three or four iterations of 
Newton's method. 


Validation 


Chien's model [5], as described, was implemented directly into the RNS3D code, borrowing 
heavily from the coding of the stream function-vorticity equations. That is, the k and e equations 
were coupled and solved using the 2x2 alternating-direction implicit (ADI) solver which was 
already in place. 

Initial solutions were obtained in a straight, circular duct at a Reynolds number of 150,000 
and compared with expected near-wall behaviors. The results are shown in Figures 1 and 2, 
attached at the end of this report. Figure 1 shows the fully-developed velocity profile for the 
circular pipe flow. It can be seen that the viscous sublayer and buffer layer are captured well, 
though the logarithmic region has a slope and intercept which may be slightly different than 
expected. Figure 2 shows the turbulent kinetic energy scaled by the friction velocity squared. 
Its near-wall behavior is again as expected and fits between the empirical bounds. In the 
logarithmic region, the turbulent kinetic energy levels off to a near-constant value; however, the 
value attained is slightly higher than what is typically expected. 

An extensive review of solutions of this type is available by Hreynya, et. al. [11]. 
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III. ANISOTROPIC TURBULENCE MODELING 


Two-equation turbulence models, especially those which depend on k and e, have become 
popular for modem engineering applications of computational fluid dynamics. They are more 
general than algebraic mixing-length models and are much less computationally expensive than 
full Reynolds stress models (i.e., those which model and solve each component of the Reynolds 
stress tensor separately). One drawback, however, is the assumption of an isotropic eddy 
viscosity, which causes the normal Reynolds stresses in most cases to be roughly, but 
erroneously, equal. 

As early as 1957, Rivlin had noted the potential importance of normal Reynolds stress 
differences [24]. Rivlin compared turbulent with non-Newtonian, visco-elastic flows in a non- 
circular straight duct where both were characterized by secondary velocities in a plane normal 
to the primary flow direction. Noting that the secondary flow in the non-Newtonian fluid was 
generated by differences in normal stress, Rivlin suggested that differences in the normal 
Reynolds stresses might explain the secondary motion in turbulent flow. Several years later, 
Gessner and Jones showed experimentally that this was indeed the case [9], Building on this 
work, researchers such as Lumley suggested that a nonlinear constitutive relation for the 
Reynolds stresses, comparable to that used for non-Newtonian fluids, might be appropriate [19]. 
These suggestions gave rise to two schools of nonlinear turbulence modeling in which researchers 
hoped to capture anisotropic effects with a nonlinear, or "anisotropic," eddy viscosity without 
resorting to the solution of six second-order, nonlinear, partial differential closure equations. 

The first class of models, the "algebraic Reynolds stress models" (ARSM), such as that 
developed by Rodi (1976), are typically derived from the Reynolds stress equations by relating 
the transport of the Reynolds stress to the transport of the turbulent kinetic energy [25,26]. The 
result is an implicit, nonlinear, algebraic relation for the Reynolds stresses and is also a function 
of k and e. 

A second class of models, largely credited to Speziale (1987), has been termed the "nonlinear 
k-e models" [30], These models, though they have been derived in different ways by different 
investigators, all assume that the relationship between the Reynolds stress and the mean strain 
rate tensor is not linear but includes terms which are quadratic in the mean strain rate tensor as 
well. This idea, which provides an explicit relationship for the Reynolds stresses, is not new, but 
instead has been applied before in visco-elastic fluid dynamics and rarefied gas dynamics 
modeling [24,39]. 

It is only in the past few years that a clear understanding of the relationship between these 
two models was achieved in light of an exact explicit solution to Rodi's implicit model [8]. The 
so-called "explicit ARSM," then, may represent the state-of-the-art in anisotropic modeling and 
helps to clearly define the relationship among the algebraic models described above. 

Implicit Algebraic Reynolds Stress Modeling 

In deriving an implicit nonlinear stress-strain relation for the Reynolds stresses, x if Rodi [25,26] 
began with the exact equations for the Reynolds stresses and the turbulent kinetic energy and 
progressed with the assumption that their transport is proportional in the following way: 
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Next, Rodi used Kolmogorov's hypothesis of local isotropy for the dissipation in which 
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The pressure strain can be broken into the sum of its slow and rapid parts if buoyancy and 
wall effects are neglected. For the slow pressure strain, Rotta's linear retum-to-isotropy 
hypothesis (1972) was invoked, while for the rapid pressure strain, a model developed by Naot, 
Shavit, and Wolfshtein (1970) or W.C. Reynolds (1970) can be used [17]. The resulting model 
for the pressure strain has 
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where 1.4 ^ C, ^ 2.2, and .60 $ y <, .78. 

The nonlinear model of Rodi is obtained by substituting Eqs. (22)-(23) into Eq. (20), which 
gives 

n ( -7 ) _ 

(24) 

Typically, this implicit model is implemented by lagging the production, P , and then inverting 
the linear system of equations for the Reynolds stress components. 

The model gives improved results for the square duct problem, where turbulence-driven 
secondary flow is induced [6,26], and has been a popular improvement in flows which contain 
swirl or rotational effects [2,15,34], 
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Nonlinear k-e Modeling 

Despite some success with the implicit ARSM, some researchers have complained that the model 
is based on "drastic assumptions which have not yet been verified" [21], while others added that 
the implicit model is "cumbersome to implement in complex flows" and that "numerical stiffness 
problems can result from the need for successive matrix inversions at each iteration" [8]. Thus, 
an explicit nonlinear model was sought. 

Assuming that the Boussinesq approximation represents only the leading terms of a series 
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expansion about the strain rate, Speziale derived a nonlinear k-e model by means of a quadratic 
expansion subject to physical constraints such as dimensional and tensorial invariance, 
realizability, and material frame-indifference [30]. Independently, Yoshizawa used a multiscale 
Direct-Interaction Approximation to develop a nonlinear k-e model which included thermal 
buoyancy effects and an anisotropic model for turbulent heat flux as well [35,36]. Several 
subsequent works refined these models (including the addition of near-wall modifications) using 
two- and three-dimensional channel flows [20,21,22,23,27], More recently, Shih and Lumley 
developed a nonlinear k-e model, again using an expansion through terms quadratic in the strain 
rate and subjecting it to constraints which make it fully realizable with coefficients which depend 
on the strain invariants [29]. 

The various constant-coefficient nonlinear k-e models bear a striking similarity which allows 
them to be represented in common form. Each model proposes 
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where P t = 2/3 p k is the turbulent pressure, S ^ represents terms which are quadratic in the mean 
velocity gradients, C p is a vector of empirical constants, and F l} is a general, second-order tensor. 
Summation is performed over p. and so S pi/ takes three forms: 
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The conventional linear k-e model, in which the Reynolds stress is modeled by the Boussinesq 
approximation, is obtained when C p and F t] are zero. The major differences in the models can 
be confined to the trailing term. F tj . 

Referring to Table 1, the 1984 model by Yoshizawa included a material derivative of the 
strain rate and third-order second derivative term, but they were never recommended for use. 
The 1986 model by Yoshizawa includes a term to model thermal buoyancy effects that was used 
to obtain good results for the flow over a heated flat plate [28], but it seems to be frequently 
overlooked by other researchers [7,10]. 

Speziale's 1987 model includes a material derivative of the strain rate which was felt to be 
significant in his calculation of the flow over a backward facing step [32], but it has more 
recently been proven to be unstable. The instability was investigated by the current authors by 
analogy to the Burnett relation for laminar stress, and the term has been dropped in practice, both 
by the current authors and by Speziale himself. 

The model by Myong and Kasagi includes a term for near-wall effects in a two-dimensional 
channel. For that term, m and n represent the streamwise and normal directions in the two- 
dimensional channel, and no summation is performed over them. As a result, the term is not 
generally applicable [20]. 

The constants in Eq. (25), listed in Table 1, are empirical. It can be seen that agreement 
among the various investigators for C, and C, is fairly good. This may be attributed to the fact 
that most models were calibrated, either experimentally or computationally, with two-dimensional 
channel flow in which C 2 does not appear. The agreement on C 2 is very poor, apparently due 
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to the relative lack of effort spent on its calibration. 

Much of the nonlinear k-e modeling effort has focused on two-dimensional problems, such 
as backward facing steps with and without heat transfer [7,16,30,32,33], heated flat plates [28], 
combustor flows [13], and rotating channel flows [37,38], Applications in three-dimensional 
problems, on the other hand, have been largely unaddressed in the literature. Speziale applied 
his model in a curved, developing duct flow and expressed optimism over the results [12]; 
however, the turbulence quantities which were used in the anisotropic eddy viscosity were 
obtained from experimental data, not from the turbulence equations. 

Explicit Algebraic Reynolds Stress Modeling 


The similarity between the implicit ARSM and explicit nonlinear k-e models can be easily 
shown. Toward that end, an approximate solution of the implicit ARSM is obtained by assuming 
that the Reynolds stress can be represented by the series expansion 

x = -P 6 +A \iE +••• (27) 

ij t ij l'r ij y ’ 


where higher-order terms (squares, cubes, etc.) in E {j and Q, r the irrotational and rotational strain 
rate tensors, are included but not shown. This series expansion is substituted for the Reynolds 
stresses appearing in Eq. (24), and terms which are higher than second order in the strain rate 
tensors are subsequently neglected. The resulting expression is 
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It is expected that A, = 2; however, C M = .09, .6 < y < .78, and 1.4 < C, < 2.2 imply A, > 2.72. On 
the other hand, for the purposes of comparison, choosing the values y = 0.8 and C t = 2.48, gives 
A, = 2, A 2 = 4, and A 3 = 2. This explicit form of Rodi's model can also be described using Table 
1, where its resemblance to the explicit models discussed above can be seen. 

Furthermore, recent research, and the availability of computer-aided symbolic manipulation, 
have allowed for the exact solution of a generalized form of the implicit ARSM, where a range 
of linear pressure strain models and non-inertial effects are included [8], 

The general solution for three-dimensional flow is quite complicated, but in the limit of two- 
dimensional flow, the solution can be written 
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where e jJk is the permutation tensor, and go, is the system rotation tensor. The term p t * uses the 
standard definition for eddy viscosity but with C M replaced by Cj which is now a function of 
the ratio of production to dissipation. (For homogeneous flows in equilibrium, Cj ~ .09.) In 
addition, each of the four B terms is a constant multiple of 
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The right-hand inequality in Eq. (32) holds for small rates of irrotational strain, or when the time 
scale of the turbulence is much smaller than the time scale set by the irrotational strain rate. 
When q 2 is also small, B 0 » 1. Then, the B terms become constant and can be set, using the 
Launder, Reece, and Rodi model constants, to be B i ~ 2, B 2 ~ 1.25, B 3 = 1.7, and B 4 ~ 6.7. 

Thus, under the assumption of homogeneous turbulence in equilibrium, in the limit of two- 
dimensional flow and for small rates of strain, the explicit ARSM of Gatski and Speziale can be 
included in Table 1 in the form of the nonlinear k-e model, where the Launder, Reece, and Rodi 
model has been used to set the constants. 

It is more important, however, to note that for sufficiently large rates of strain, the term in 
Eq. (31) can become singular. In this explicit form, the term can be regularized so that it is well- 
behaved in these situations where the model is not formally valid; however, Eq. (31) highlights 
why the implicit ARSM has been found, at times, to be unstable [8]. 


Discussion 


In summary, while Rodi's implicit ARSM was a major step in anisotropic modeling, the explicit 
solution of Gatski and Speziale is more convenient and allows for the removal of any inherent 
singularities. In the limit of two-dimensional flow, the explicit solution reduces to a quadratic 
form which can represented by any of the nonlinear k-e models in the limit of small rates of 
strain. For many flows, these limits are not severe, and so the nonlinear k-e models will capture 
most of the effects available from the original ARSM. 

The current authors, for instance, have obtained improved results for internal duct flows using 
a constant-coefficient nonlinear k-e model. The average fluctuating velocities in fully-developed 
two-dimensional channel flow, which would be predicted equal by a linear model, are shown in 
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Figure 3 compared to experimental data, where it can be seen that the nonlinear model provides 
considerable improvement. Additionally, the secondary velocities in a fully-developed three- 
dimensional square duct flow were obtained and are shown in Figure 4. Here, the comparison 
to data shows that the velocities are much improved in contrast to the zero secondary-velocity 
that the linear model must give. And while highly three-dimensional flows have remained largely 
unaddressed in the literature, the current authors have obtained improved normal Reynolds 
stresses for the flow over a 6:1 prolate spheroid at 10° angle-of-attack, the results of which have 
been submitted for publication. These results suggest that the nonlinear k-e model gives 
improvements even in the case of highly three-dimensional flow; however, knowing that 
nonlinear k-e models represent ARSMs in the limit of two-dimensional flow, such applications 
should be regarded with considerable caution. 

For flows of interest to RNS3D, the nonlinear k-e model should offer many advantages. Any 
of the generic constant-coefficient models should suffice, but particular attention should be paid 
to the work of Myong and Kasagi [20.2 1 ] as it seems to be the only effort to study anisotropic 
effects in the near-wall region, a region of consequence for the current effort. As mentioned, the 
near-wall term added by Myong and Kasagi was not frame invariant, but may be suitable under 
the RNS3D coordinate system or amenable to generalization. 

Finally, the model of Shih and Lumley [29], which had not been studied by the current 
authors until recently, appears to be a viable alternative to the models discussed in more detail 
here. Clearly, in satisfying realizability, the model becomes slightly more complicated. In any 
case, the near-wall region would still need separate attention as it has not been addressed by Shih 
and Lumley [29]. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


Conclusions 

The current effort has incorporated a two-equation k-e model with low-Reynolds number 
modifications in order to address deficiencies in RNS3D in the prediction of duct flows with 
vortex generators. Coding and implementation is entirely consistent with RNS3D coding and 
assumptions, which will allow NASA LeRC personnel already familiar with the code some 
flexibility with regards to modifications. The equations appear to be efficient in solving a 
circular pipe flow and capture the near-wall region well. 

Proposed Future Activity 

The goal of improving the capability of NASA LeRC personnel in predicting the turbulence 
quantities in duct flows with vortex generators has been largely realized in the current effort. 
The RNS3D code now includes a more general two-equation turbulence model applicable in the 
near-wall region, and it exists in a form which will be familiar to RNS3D users. Future efforts 
should include careful examination of the coding to additional applications to ensure that the 
model is made suitably robust. The ability of NASA LeRC personnel to use the new model will 
greatly expedite that process. 

Nonlinear k-e models, such as those proposed by Yoshizawa (1984) [35] and Speziale (1987) 
[30], have become an attractive option to users of the k-e model, and more recent work by Gatski 
and Speziale has shown the extent to which these constant coefficient, explicit nonlinear models 
anticipated the solution of the implicit algebraic models such as that by Rodi [8]. But it has also 
shown that the models only formally apply in the limit of two-dimensional flow. In addition, 
while most models do not take the near-wall region into account, work by Myong and Kasagi 
[20] includes a near-wall term which, with slight modification, may be suitable for use in duct 
flows controlled by vortex generators. Because the RNS3D code now includes a two-equation 
model, the step to include anisotropic effects becomes a relatively minor one, and its 
implementation should be included in any future activity. 
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Table 1: A Comparison of Nonlinear k-e Models 
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